Emergence of glass-like dynamics for dissipative and strongly interacting bosons 



o 

(N 
O 

Q 



Dario Poletti, 1 Peter Barmettlcr, 2 Antoinc Georges, 3 ' 4, 5 and Corinna Kollath 2 

Singapore University of Technology and Design, 20 Dover Drive, 138682 Singapore 
2 Departement de Physigue Theorique, Universite de Geneve, CH-1211 Geneve, Switzerland. 

College de France, 11 place Marcelin Berthelot, 75005 Paris, France. 
"'Centre de Physique Theorique, Ecole Poly technique, CNRS, 91128 Palaiseau Cedex, France. 
5 DPMC-MaNEP, Universite de Geneve, CH-1211 Geneve, Switzerland. 

We study the dynamics of a strongly interacting bosonic quantum gas in an optical lattice potential 
under the effect of a dissipative environment. We show that the interplay between the dissipative 
process and the Hamiltonian evolution leads to an unconventional dynamical behavior of observables, 
such as local number fluctuations. In particular we show, both analytically and numerically, the 
emergence of an anomalous diffusive evolution at short times and, at long times, a glass- like dynamics 
dominated by rare events. This complex dynamics reveals information on the level structure of the 
strongly interacting gas. 
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Unconventional, non-exponential, relaxation dynam- 
ics of a perturbed system towards equilibrium has at- 
tracted a lot of interest over decades. Already in 1847, 
Kohlrausch [l[ observed a stretched exponential decay 
in time t, i.e. e -1 -'/ 1 " ' with a G (0, 1) and to a posi- 
tive constant, of the discharge of capacitors fabricated 
from glasses. Since then, such a decay has been observed 
in many systems such as molecules and polymersf^, Q , 
spin glasses 0, Q nanosized magnetic particles 0| and 
certainly amorphous silicon 0, H| . 

A broad variety of theoretical approaches has been de- 
veloped to explain the mechanism of unconventional re- 
laxation dynamics [1- 11|. In many of these approaches, 
e.g. the treatment of the Griffiths phase in disordered 
spin systems [l2[ , rare configurations have been identified 
to play a key role. These configurations have an expo- 
nentially small probability to occurr, and therefore con- 
tribute minimally to the short-time dynamics. However, 
because their relaxation time scale is very long, these 
rare configurations may dominate the long-time evolu- 
tion. This can result in subexponential relaxation, often 
a stretched exponential decay, which we refer to as 'glass- 
like' in the following. 

In this work, we uncover that also in seemingly sim- 
ple and quantum many body systems without disorder, 
as the Bose-Hubbard model, the dissipative coupling to 
a Markovian, i.e. memory less, environment can cause 
glass-like dynamics. Wc show that heating up such a 
clean system by an external coupling can force rare con- 
figurations to be visited. These are the configurations for 
which one site has a relative large atom number. Since 
the transition to these rare configurations arc associated 
to long time-scales due to the presence of strong interac- 
tions, they are at the origin of an unconventional dynam- 
ics of stretched exponential form. Thus we propose that 
the heating dynamics of a non-random quantum many- 
body system in the presence of dissipation can display 
glassy behavior, which reveals the complex structure of 



its configuration space and energy spectrum. 

A non-exponential decay is in contrast to the typical 
evolution found for quantum many body systems cou- 
pled to an environment. In these many body systems, 
the decay is often dominated by an exponential dynam- 
ics, as e.g., the counterintuitive Zeno effect or the 
relaxation to a desirable state driven by an artificially en- 



gineered environment [171 . [18[ . Only recently, first signs 



of an intriguing slowing down of the heatin g dy namics for 
interacting bosonic 13, 2(| and fermionic }2ll ] gases and 
an algebraic decay for a specially designed environment 
which imprints coherence [22 ] have been predicted. How- 
ever, up to now, the understanding of the large variety of 
dynamical behaviors that occur in an interacting many 
body system and of its origin is still a great challenge. 

In this work, we study the heating of ultracold bosonic 
atoms in an optical lattice. The unitary dynamics of the 
bosons can be well described by the Bose-Hubbard model 
23|, 24 1 . The heating can, for example, be due to spon- 
taneous emissions of the red detuned optical lattice light 
field itself which has been identified as one of 

the most important sources of heating in optical lattices. 
We show that over a large range of a rescaled time r, 
the heating process can be well described by a diffusion 
equation of the form 



dp{x, ■ 
dr 



d_ 

dx 



D(x,t 



dp(x, t) 
dx 



■ F(x,t)p(x,t) 



(1) 



where p(x,r) is the probability that one lattice site has 
an atom number occupation xf, where / is the aver- 
age filling of the lattice. Interestingly, in the present 
context, the diffusion function D(x,t) and force field 
F(x, t) depend on the distribution function p(x, r) it- 
self. This dependence makes the equation above a non- 
linear diffusion equation in which D and F have to be 
determined self-consistently. Fortunately, the behavior 
of D and F at both the short-times and long-times can 
be obtained analytically. At short times, we will show 
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that the diffusion function is dominated by a divergence 
D oc l/(x — l) 2 close to i « 1 while the force is neg- 
ligible thus, causing an anomalous diffusive process. In 
contrast, for large times, and for x large enough, the dif- 
fusion function decreases as D{x, r) sa 1/x and the force 
term F(x,t) ps corresponds to an effective loga- 

rithmic potential. Here, the occupation of rare states 
at large x dominates the long time evolution. As a re- 
sult, the relaxation to the steady state slows down and 
the fluctuations k = (nj) — (fij) 2 obey a stretched ex- 
ponential decay as shown in Fig. [T] This glass-like dy- 
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FIG. 1: (color online) 1 — k(t)//c(oo) vs square-root of rescaled 
time t for U/fry = 10. Numerical results of eq. (|3]) are shown 
for different fillings, / = 0.1,0.5, 1,3 in the direction of the 
arrow, as solid lines and corresponding stretched exponential 
fits as dashed lines. The analytical result (0 of the diffusion 
equation (JXJ) is shown as a black dot-dashed line. 

namics is typical for disordered systems, however here it 
arises, during the heating process, purely from the inter- 
play of kinetic and interaction energy. Via the diffusion 
equation ([TJ, unexpected connections to other intriguing 
physical problems arise. One example is the emergence 
of non-ergodic and super-aging behavior for diffusion in 
a logarithmic potential, but with a constant D (2(| . This 
situation is realized in another category of experiments 
with dissipative cold atoms in optical lattices j27l |28|. 

The dissipative system under study, composed of N 
bosonic atoms in a deep optical lattice of N s sites with 
filling / = N/N s and connectivity (number of nearest 
neighbors per site) z, can be described by the following 
master equation [1S|, |25j : 



d t p = --[H,p] 



V(p) 



(2) 



The first term describes the unitary evolution of the 
density matrix p. This evolution is governed, in the 
single Bloch band limit, by the Hamiltonian H = 
~^E(j,() + f Ej % ("i _ !) where (J, I) denotes 
pairs of neighboring sites. The operators bj and bj arc 
bosonic creation and annihilation operators on site j and 
fij = bjbj counts the number of atoms. The second term 



models the dissipative coupling to a Markovian environ- 
ment, as the optical lattice itself 25 1, via the local 



density with strength 7. We have restricted the descrip- 
tion to the lowest Bloch band of the optical lattice poten- 
tial. The validity of this approximation will be discussed 
in the conclusions. 

In the following, we study in detail the heating dy- 
namics of a system initially in its ground state state with 
respect to H under the joint action of dissipation and 
H. We concentrate on the strongly interacting regime 
U ^> J,hy. The dissipator causes the off-diagonal ele- 
ments of the density matrix, in the Fock basis, to decay 
towards the decoherence free subspace which consists of 
all possible diagonal density matrices p. In the presence 
of the hopping term, the heating process drives the sys- 
tem to a unique steady state p(t = 00) = -hi, the high- 
est entropy state [2(J. Here, M is the dimension of the 
Hilbert space at fixed atom number N, and I the iden- 
tity operator. The approach of this infinite temperature 
state, can be described for jt 3> 1, by adiabatically elim- 
inating [29| the small off-diagonal elements. A closed set 
of classical rat e eq uations for the diagonal elements of p 
is obtained 2(], 30| . The diagonal configurations that are 
connected are those for which a particle is moved from 
a site with occupation ml to one of its neighbors with 
occupation to. The process occurs via virtual hopping 
to and from an off-diagonal element of the density ma- 
trix. The corresponding transition rates emerge from the 
imaginary part of the virtual hopping transition coeffi- 
cients, i.e. Im\-i — tP/yft-, ■ fc ) = ttT(to,to',1) with 

f 2 

T(m, to', d) = 



{m-\-X—m l )U-\-ih'y 
(m+6 d<1 )(m'+8 d 



1) 



(m-m'+d) 2 + (fi, 7 /C/) 2 



and t* = To 



study this dynamics, we use a separable and translation- 
ally invariant ansatz p(t) = \%2 n p(n, t)\n) (n\] where 
j runs over all the lattice sites and n over all the possi- 
ble occupations of each site. The probability distribution 
p(n, t) of the single site occupation evolves as 



9 T p(n,r)= T{m,n,d) 



m,d—±l 



p(m — d, T)p(n + d, r) 



p(m,r)p(n,T) 



(3) 



of equation j2]), V (p) = 7^, 



ijprij 



\h)p - 



\phfj 



where r = t/t*. In the following, we will use this equation 
to obtain numerically the evolution of the system. An 
impression of a typical evolution of the occupation num- 
ber distribution can be acquired by studying Fig. [2] At 
short times, but "ft > 1, the very narrow initial distribu- 
tion around the average filling / broadens symmetrically 
(see inset of Fig. [2]) . After the rapid broadening a new 
regime with an asymmetric evolution sets in, in which the 
tail of the distribution slowly converges towards the ex- 

1 ( f \ n+1 

pected asymptotic distribution p(n, 00) = j ( y+/J 
p(n, 00) is the single site reduced density matrix of the 
full asymptotic density matrix p(t = 00) = -jjX (3l| . 
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FIG. 2: (color online) Evolution of the density matrix ele- 
ments p(n, t) (solid lines) in a semi-logarithmic plot, versus n 
for large rescaled times r between 0.1 and 50 (not-equidistant) 
in the direction of the red arrow. The red dashed line shows 
the asymptotic limit. The inset shows the same evolution at 
shorter times r between 0.0002 and 0.1 (not equidistant) in a 
linear plot. Parameters: / = 3, U /fry = 10. 



In order to obtain analytical insight into these very 
different regimes of the evolution, we take the continuum 
limit of eq. ^ for large /. The continuous on-site occu- 
pation number distribution p(x = nj f) — fp(n) fulfills 
the non-linear integro-diffcrential equation ([TJ with 



D(x,t) 



xy p{y,r) 
(x - y) 2 + e 2 



F(x,r) = f 
Jo 



x y d v p(y, T ) 

(x — y) 2 + £' 



dy, 



dy 



(4a) 
(4b) 



and e = fvy/fU. The complicated form of D(x,t) and 
F(x, t) stems from the configuration dependent rates 
and triggers a wide range of rich phenomena. Note 
that the structure of ((4]) insures that both the total 
probability (J Q p{x 1 r)dx = 1) and the average popu- 
lation (J^° p(x,r)dx = 1) are conserved. Further, it 
can be checked that the asymptotic solution of ([1} is 
p(x, oo) = e~ x , which is the continuum limit of the steady 
state in the large-/ limit. The continuum description is 
justified for / large and finite e, assuming that p(x,r) 
varies smoothly enough on scales of the order of 1//. In 
the present case, the strongest variations of distributions 
arc due to the initial state, especially for a low filling 
/, which smooth out rather quickly and the continuum 
description is highly accurate over a wide time range. 

In the following we will solve analytically the diffusion 
equation JTJ) in the short time and in the long time limits 
focusing on the evolution of the particle distribution and 
the local density fluctuations k. 

Short time relaxation: Within the diffusion equation 
([1]), initially the distribution p is strongly peaked and 
symmetric around the value x = 1. For such a distri- 
bution the force is ncglcgiblc compared to the diffusion 



function. The diffusion equation at x « 1 can be approx- 
imated by d T p{x, t) = d x 



_ ( z i)a +e i> flxP(g, r)J . 
This leads to a dynamics which in the analytically solv- 
able limit e —> is given by an anomalous diffusion 

of the form p(x,r) = 4r(5/ 4 ) T i/ 4 e . T(s) is the 

gamma function [32| . Using this analytical solution for 
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FIG. 3: (color online) Local particle fluctuations n/ f 2 ver- 
sus rescaled time r for U/Tvy = 10. Numerical results (solid 
lines) for various filling, / = 3, 5, 7, 9, 20 in the direction of 
the arrow, are obtained solving eq. ([3]). The approximate an- 
alytical solution k/ f 2 — p(5/4) \pr of the diffusion equation is 
represented by the dashed line. 



the particle distribution, the local number fluctuations 
k(t) I f 2 = J Q (x 2 — l)p(x,r)dx exhibit a power-law re- 
laxation with n/f 2 = p(5/4) V^- This analytical result 
is in excellent agreement with the numerical results in 
Fig. [3] obtained by solving eq. Deviations arc found 
at small fillings and short times, where the approxima- 
tion e — > is not justified. However, for larger fillings, for 
example / = 20, the time-regime in which the power-law 
appears is already large. 

Physically, this very rapid initial broadening of the par- 
ticle distribution translates into the fast creation of small 
particle fluctuations around the average value caused by 
the heating. These fluctuations arise via virtual exci- 
tations of low energetic cost of order O(U) which thus 
can be reached rapidly. This behavior is similar to the 
dynamics observed in a double well potential (2p| . 

Long time relaxation: The obtained short time solu- 
tion breaks down as the distribution approaches the re- 
flective boundary at x = [33| : the distribution is no 
longer symmetric around x = 1, and the combined ac- 
tion of the force term with the diffusion will drive the 
system towards its large time asymptotics p(x, oo) = e~ x . 
Physically, the exponential suppression of large values of 
x corresponds to the rareness of the states with a high 
number of particles accumulated on a single site. 

Considering the diffusion function and the force cor- 
responding to these large values of x, one finds, in our 
study, that the dominating contributions arc indepen- 
dent of p and given by D(x,t) sa —F{x,t) -. The 
underlying quantum mechanical process behind this slow 
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diffusion is the large energy cost of the virtual states via 
which the diffusion processes at high x takes place. In 
order to describe the approach towards the asymptotic 
solution (see Fig. [2]), the ansatz p(x, r) = p(x, oo)g(x, r) 
is suitable. The evolution of the function g is shown in 
Fig. S] and suggests a scaling form g{x 7 r) = g(r)) with 
Here a and b are some functions of r to be 



'/ 



b(T 



determined. Adopting this ansatz and neglecting unim- 
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FIG. 4: (color online) (a) The blue-dotted line shows the evo- 
lution of g(x, t) versus x for large rescaled times r £ [28, 710] 
in the direction of the red arrow. Inset: Plot of g(x,r) versus 
g = (x — o(r))/6(r) (blue line) and fit with an error function 
(red dashed line), (b) a(r) is plotted versus \fr (blue line) 
and compared to a linear fit (red dashed line), (c) b(r) is 
plotted versus r 1//4 (blue line) and compared to a linear fit 
(red dashed line). Parameters: / = 0.5, U/tvy = 15. 



portant terms, equation (jTJ) reduces to 

dg d 2 g 



b (da — !) + ?/( abb + b"d 



dg dg 2 



= 



(5) 



where the dot stands for the derivative with respect to r. 
This scaling equation must be independent of r for the 
scaling ansatz to be a solution. This leads to the only 
physically meaningful solution given by 



p(x,oo) J 
P[ x , ) = o < 1 - erf 



V3 x - 
T (2r)V4 



(6) 



where erf is the error function |32| . 

Figs. BJb-c) show that, at long enough times, the nu- 
merical data and the proposed analytical r-dependence 
a(r) oc yfr and b[f) oc r 1 / 4 match accurately. However, 
let us note that the numerical results still show devia- 
tions from the exact analytical prefactors. We verified 
that these deviations become smaller with increasing fill- 
ing; this leads us to conclude that the analytical form 
will be applicable in the large / limit. 



For large enough times, the obtained solution gives a 
very good approximation in the entire range of x. The 
reason for this are the fast relaxation times at low values 
of x. This leads to a fast initial relaxation and, there- 
after, only small relative changes occur at small x. These 
changes in the probability distribution are mainly con- 
nected to the variations at large x via particle number 
conservation. Thus, we can use the obtained solution 
to calculate the local particle fluctuations K. The corre- 
sponding integral can be solved analytically giving 



k(t) oc h(r) e 



(7) 



where h depends algebraically on r. We thus have shown 
analytically the emergence of the stretched exponential 
behavior. This finding, as depicted in Fig. [TJ compares 
well to the numerical solution of eq. ((3|). However, the 
time at which the stretched exponential occurs, and the 
detailed decay, depend on the filling. In particular, the 
stretched exponential occurs later for larger fillings. 

The origin of the stretched exponential evolution can 
be traced back to the occurrence of rare events in the 
large — x tail of the distribution function e~ x which are 
connected to decreasing time scales oc 1/x. This is clearly 
seen using a saddle point argument for the evaluation of 
the fluctuations employing only these two ingredients, i.e. 



f 



oo 
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The saddle point is determined by 
and leads to the main contribution 



Bose- 



K 

dx [ X T ~x) 1*0 

ft w e _v/ ^* recovering the stretched exponential. 

To summarize, we have uncovered, in the 
Hubbard model coupled to a dissipative environment, 
two different regimes of unconventional relaxation: at 
short times and large enough fillings, a power-law regime, 
while at large times, and any filling, a stretched exponen- 
tial regime. This last regime is dominated by rare events 
which correspond to the occupation of a single site with a 
large number of atoms. Thus, very slow transition rates, 
due to the high energetic cost of the processes connecting 
these rare configurations, dominate the long time dynam- 
ics. This emergent glass-like dynamics is thus a signature 
of the level spacing between energy bands of the Bose- 
Hubbard Hamiltonian and of the number of configura- 
tions present in each band. We would like to point out 
that this realization of a glass regime is very natural and 
simple: the origin of the Kohlrausch relaxation is not due 
to disorder or to artificially imposed constraints |34f but 
to the structure of the energy levels of the system itself. 

An experimental observation of these relaxation 
regimes, especially of the glass-like relaxation in a non- 
random system, would be fascinating. This requires an 
optimisation for the choice of the filling: for larger fill- 
ings the algebraic decay is more pronounced, whereas 
for a low filling the glass-like behavior sets in already at 
shorter time-scales. Additionally, one has to assure that 
the influence of higher bands is not overwhelming the 
effect. Particular care has to be taken that the energy 
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of the created excitations, e.g. the configuration with a 
large number of atoms accumulated at a certain lattice 
site, is still lying well below the next Bloch band. 
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